setwd("~/Desktop/Research/NSF/WCT hybridization/WCT_RBT_Genomics/Long_results")
RBT_frequencies <- read.delim("~/Desktop/Research/NSF/WCT hybridization/WCT_RBT_Genomics/Long_results/ConGen_LongTest_data.txt")
data <- RBT_frequencies
L <- length(data$FinRBT)  #Number of loci
M <- mean(data$FinRBT)    #Average admixture across loci
###perform calculations
Vinv <- matrix(0,L,L)
diag(Vinv) <- 1/(M*(1 - M))
expvar <- mean(data$FinRBT*(1 - data$FinRBT))
MSE <- t(data$FinRBT - M)%*%Vinv%*%(data$FinRBT - M)/(L - 1)
MSEvec <- MSE + expvar*(1/(2*data$FinN) - 1/(2*mean(data$FinN)))
VMi <- data$FinRBT*(1 - data$FinRBT)*MSEvec
###obtain statistical results by locus
data <- transform(data, Fin_chi = (data$FinRBT - M)^2/VMi) #chi-square locus component
data <- transform(data, Fin_resid = (data$FinRBT - M)/sqrt(VMi)) #residual chi-square
data <- transform(data, Fin_effect = (data$Fin_resid)/sqrt(data$FinN)) #effect size
data <- transform(data, Fin_p = pchisq(data$Fin_chi,1,lower.tail=FALSE))  #p-value by locus
###write results to data table
write.table(format(data,scientific=FALSE),file="long_output.txt",quote=FALSE,row.names=FALSE,sep="\t")

L <- length(data$CycCRBT)  #Number of loci
M <- mean(data$CycCRBT)    #Average admixture across loci
###perform calculations
Vinv <- matrix(0,L,L)
diag(Vinv) <- 1/(M*(1 - M))
expvar <- mean(data$CycCRBT*(1 - data$CycCRBT))
MSE <- t(data$CycCRBT - M)%*%Vinv%*%(data$CycCRBT - M)/(L - 1)
MSEvec <- MSE + expvar*(1/(2*data$CycCN) - 1/(2*mean(data$CycCN)))
VMi <- data$CycCRBT*(1 - data$CycCRBT)*MSEvec
###obtain statistical results by locus
data <- transform(data, CycC_chi = (data$CycCRBT - M)^2/VMi) #chi-square locus component
data <- transform(data, CycC_resid = (data$CycCRBT - M)/sqrt(VMi)) #residual chi-square
data <- transform(data, CycC_effect = (data$CycC_resid)/sqrt(data$CycCN)) #effect size
data <- transform(data, CycC_p = pchisq(data$CycC_chi,1,lower.tail=FALSE))  #p-value by locus
###write results to data table
write.table(format(data,scientific=FALSE),file="long_output.txt",quote=FALSE,row.names=FALSE,sep="\t")

L <- length(data$HayRBT)  #Number of loci
M <- mean(data$HayRBT)    #Average admixture across loci
###perform calculations
Vinv <- matrix(0,L,L)
diag(Vinv) <- 1/(M*(1 - M))
expvar <- mean(data$HayRBT*(1 - data$HayRBT))
MSE <- t(data$HayRBT - M)%*%Vinv%*%(data$HayRBT - M)/(L - 1)
MSEvec <- MSE + expvar*(1/(2*data$HayN) - 1/(2*mean(data$HayN)))
VMi <- data$HayRBT*(1 - data$HayRBT)*MSEvec
###obtain statistical results by locus
data <- transform(data, Hay_chi = (data$HayRBT - M)^2/VMi) #chi-square locus component
data <- transform(data, Hay_resid = (data$HayRBT - M)/sqrt(VMi)) #residual chi-square
data <- transform(data, Hay_effect = (data$Hay_resid)/sqrt(data$HayN)) #effect size
data <- transform(data, Hay_p = pchisq(data$Hay_chi,1,lower.tail=FALSE))  #p-value by locus
###write results to data table
write.table(format(data,scientific=FALSE),file="long_output.txt",quote=FALSE,row.names=FALSE,sep="\t")

##Calculate combined P value using Fisher's combined probability test

Long_results <- read.delim("~/Desktop/Research/NSF/WCT hybridization/WCT_RBT_Genomics/Long_results/long_output.txt")
data <- Long_results
n <- 3 ##number of P values being combined
df <- 2*n
data <- transform(data, Combined_P = pchisq( -2*(log(data$Fin_p)+log(data$CycC_p)+log(data$Hay_p)), df, lower.tail=FALSE))
write.table(format(data,scientific=FALSE),file="long_output.txt",quote=FALSE,row.names=FALSE,sep="\t")

###Correction for multiple tests
Combined_P <- pchisq( -2*(log(data$Fin_p)+log(data$CycC_p)+log(data$Hay_p)
fdr <- p.adjust(Combined_P, "BH")
sum(fdr < 0.05)
sum(Combined_P < 0.05)
sum(Combined_P < 0.05/9380)